```{r}

library(readxl)
library(dplyr)
library(ggplot2)
library(stargazer)
library(ggpubr)
library(pwrss)
library(kableExtra)
library(knitr)



```

```{r}
############# Reading Data #################

survey = read_excel("survey_JoP.xlsx", sheet = "Sheet1")

```

```{r}
########## Figure 1: Earmark Distribution

# Discarding no response
pork = filter(survey, Pork != "N/R")

# Recoding the answers & adding new column - here we we create a new variable that corresponds to the legislator preference on an issue-area
pork <- pork %>%
  mutate(Pork_num = recode(Pork,
                             "Em certas circunstâncias, o STF tem o poder de decidir sobre a distribuição das emendas parlamentares e de revogar atos do Poder Executivo." = "Active Court", 
                             "O STF NÃO tem o poder de decidir sobre a distribuição das emendas parlamentares e NEM de revogar atos do Poder Executivo." = "Passive Court"))

# Calculating the Confidence Interval 

## Creating a count and percentage table
tabela_cenarios = pork %>% group_by(Pork_num) %>% summarise(count = n())
perct = tabela_cenarios$count/nrow(pork)
tabela_cenarios = cbind(tabela_cenarios, perct)

# Percentages - here we create objects [c1 and c2] that contains the preference percentages 
c1 = tabela_cenarios[1,3]
c2 = tabela_cenarios[2,3]

# Confidence intervals - here we calculate lower/upper confidence intervals
c1+c2
se = sqrt((c1*c2)/nrow(pork))
CI = 1.96 * se

lower.CI = c(c1-CI, c2-CI) 
upper.CI = c(c1+CI, c2+CI) 

# Confidence interval table - here we build the confidence interval table
tabela_cenarios_pork = cbind(tabela_cenarios, lower.CI, upper.CI)
tabela_cenarios_pork

# Printing the table
stargazer(tabela_cenarios_pork,summary = F, type = "text", digits=4, no.space=T, flip=F)

# Saving confidence interval table [Earmark Distribution]
writexl::write_xlsx(tabela_cenarios_pork, "Earmark_Distribution.xlsx")

#### Figure 1: Earmark Distribution

a = read_xlsx("Earmark_Distribution.xlsx")

graph_pork = ggplot(a) +
  geom_bar( aes(Pork_num , y=perct), stat="identity", fill="white", width=0.4, colour = "black", alpha= 2.5) +
  geom_errorbar( aes(x=Pork_num, ymin=lower.CI, ymax=upper.CI), width=0.1, colour="black", alpha=0.9, size=1) + ggtitle("Figure 1: Earmark Distribution") +
  xlab("") +
  ylab("") +
  theme(plot.title = element_text(family = "serif", size = 20, face = "bold", hjust = 0.5), axis.title.x = element_text(family = "serif"), axis.text = element_text(family = "serif", colour="black", size = 19) ) 

ggsave("Figure 1 - Earmark Distribution.png", width = 2105 , height = 2105, units="px", dpi= 420)

```

```{r}
########## Figure 2: Social Values

# Discarding no response
check = filter(survey, Check != "N/R")

# Recoding the answers & adding new column - here we we create a new variable that corresponds to the legislator preference on an issue-area
check <- check %>%
  mutate(Check_num = recode(Check,
                             "Em certas circunstâncias, o STF tem o poder de decidir sobre temas de comportamento (aborto e casamento) e de revogar atos do Poder Executivo." = "Active Court",
                             "O STF NÃO tem o poder de decidir sobre temas de comportamento (aborto e casamento) e NEM de revogar atos do Poder Executivo." = "Passive Court"))


# Calculating the Confidence Interval 

## Creating a count and percentage table
tabela_cenarios = check %>% group_by(Check_num) %>% summarise(count = n())
perct = tabela_cenarios$count/nrow(check)
tabela_cenarios = cbind(tabela_cenarios, perct)

# Percentages - here we create objects [c1 and c2] that contains the preference percentages 
c1 = tabela_cenarios[1,3]
c2 = tabela_cenarios[2,3]

# Confidence intervals - here we calculate lower/upper confidence intervals
c1+c2
se = sqrt((c1*c2)/nrow(check))
CI = 1.96 * se

lower.CI = c(c1-CI, c2-CI) 
upper.CI = c(c1+CI, c2+CI) 

# Confidence interval table - here we build the confidence interval table
tabela_cenarios_check = cbind(tabela_cenarios, lower.CI, upper.CI)
tabela_cenarios_check

# Printing the table
stargazer(tabela_cenarios_check,summary = F, type = "text", title="Social Values", digits=4, no.space=T, flip=F)

# Saving confidence interval table [Social Values]
writexl::write_xlsx(tabela_cenarios_check, "Social_Values.xlsx")

#### Figure 2: Social Values

a = read_xlsx("Social_Values.xlsx")

graph_check = ggplot(a) +
  geom_bar( aes(Check_num , y=perct), stat="identity", fill="white", width=0.4, colour = "black", alpha= 2.5) +
  geom_errorbar( aes(x=Check_num, ymin=lower.CI, ymax=upper.CI), width=0.1, colour="black", alpha=0.9, size=1) + ggtitle("Figure 2: Social Values") +
  xlab("") +
  ylab("") +
  theme(plot.title = element_text(family = "serif", size = 20, face = "bold", hjust = 0.5), axis.title.x = element_text(family = "serif"), axis.text = element_text(family = "serif", colour="black", size = 19) ) 

ggsave("Figure 2 - Social Values.png", width = 2105 , height = 2105, units="px", dpi= 420)

```

```{r}
########## Figure 3: Electoral Legislation

# Discarding no response
elec = filter(survey, Election != "N/R")

# Recoding the answers & adding new column - here we we create a new variable that corresponds to the legislator preference on an issue-area.  
elec <- elec %>%
  mutate(Election_num = recode(Election,
                             "Em certas circunstâncias, o STF tem o poder de decidir sobre legislação eleitoral e de revogar atos do Poder Executivo." = "Active Court", 
                             "O STF NÃO tem o poder de decidir sobre legislação eleitoral e NEM de revogar atos do Poder Executivo." = "Passive Court"))

# Calculating the Confidence Interval 

## Creating a count and percentage table
tabela_cenarios = elec %>% group_by(Election_num) %>% summarise(count = n())
perct = tabela_cenarios$count/nrow(elec)
tabela_cenarios = cbind(tabela_cenarios, perct)

# Confidence intervals - here we calculate lower/upper confidence intervals
c1 = tabela_cenarios[1,3]
c2 = tabela_cenarios[2,3]

# Confidence intervals - here we calculate lower/upper confidence intervals
c1+c2
se = sqrt((c1*c2)/nrow(pork))
CI = 1.96 * se

lower.CI = c(c1-CI, c2-CI) 
upper.CI = c(c1+CI, c2+CI) 

# Confidence interval table - here we build the confidence interval table
tabela_cenarios_elec = cbind(tabela_cenarios, lower.CI, upper.CI)
tabela_cenarios_elec

# Printing the table
stargazer(tabela_cenarios_elec,summary = F, type = "text", digits=4, no.space=T, flip=F)

# Saving confidence interval table [Electoral Legislatoin]
writexl::write_xlsx(tabela_cenarios_elec, "Electoral_Legislation.xlsx")

#### Figure 3: Electoral Legislation

a = read_xlsx("Electoral_Legislation.xlsx")

graph_elec = ggplot(a) +
  geom_bar( aes(Election_num , y=perct), stat="identity", fill="white", width=0.4, colour = "black", alpha= 2.5) +
  geom_errorbar( aes(x=Election_num, ymin=lower.CI, ymax=upper.CI), width=0.1, colour="black", alpha=0.9, size=1) + ggtitle("Figure 3: Electoral Legislation") +
  xlab("") +
  ylab("") +
  theme(plot.title = element_text(family = "serif", size = 20, face = "bold", hjust = 0.5), axis.title.x = element_text(family = "serif"), axis.text = element_text(family = "serif", colour="black", size = 19) ) 

ggsave("Figure 3 - Electoral Legislation.png", width = 2105 , height = 2105, units="px", dpi= 420)


```

```{r}
########## Figure 4: Privileged Jurisdiction

# Discarding no response
venue = filter(survey, Venue != "N/R")

# Recoding the answers & adding new column - here we we create a new variable that corresponds to the legislator preference on an issue-area.  
venue <- venue %>%
  mutate(Venue_num = recode(Venue,
                             "Em certas circunstâncias,  o STF tem o poder de decidir sobre o foro privilegiado e de revogar atos do Poder Executivo." = "Active Court", 
                             "O STF NÃO tem o poder de decidir sobre o foro privilegiado e NEM de revogar atos do Poder Executivo." = "Passive Court"))

# Calculating the Confidence Interval 

## Creating a count and percentage table
tabela_cenarios = venue %>% group_by(Venue_num) %>% summarise(count = n())
perct = tabela_cenarios$count/nrow(venue)
tabela_cenarios = cbind(tabela_cenarios, perct)

c1 = tabela_cenarios[1,3]
c2 = tabela_cenarios[2,3]

# Confidence intervals - here we calculate lower/upper confidence intervals
c1+c2
se = sqrt((c1*c2)/nrow(venue))
CI = 1.96 * se


# Confidence interval table - here we build the confidence interval table
lower.CI = c(c1-CI, c2-CI) 
upper.CI = c(c1+CI, c2+CI) 

tabela_cenarios_venue = cbind(tabela_cenarios, lower.CI, upper.CI)
tabela_cenarios_venue

# Printing the table
stargazer(tabela_cenarios_venue,summary = F, type = "text", digits=4, no.space=T, flip=F)

# Saving confidence interval table [Privileged Jurisdiction]
writexl::write_xlsx(tabela_cenarios_venue, "Privileged_Jurisdiction.xlsx")

#### Figure 4: Privileged Jurisdiction

a = read_xlsx("Privileged_Jurisdiction.xlsx")

graph_venue = ggplot(a) +
  geom_bar( aes(Venue_num , y=perct), stat="identity", fill="white", width=0.4, colour = "black", alpha= 2.5) +
  geom_errorbar( aes(x=Venue_num, ymin=lower.CI, ymax=upper.CI), width=0.1, colour="black", alpha=0.9, size=1) + ggtitle("Figure 4: Privileged Jurisdiction") +
  xlab("") +
  ylab("") +
  theme(plot.title = element_text(family = "serif", size = 19, face = "bold", hjust = 0.5), axis.title.x = element_text(family = "serif"), axis.text = element_text(family = "serif", colour="black", size = 19) ) 

ggsave("Figure 4 - Privileged Jurisdiction.png", width = 2105 , height = 2105, units="px", dpi= 420)
```

```{r}
########## Figure 5: Bureaucracy

# Discarding no response
tech = filter(survey, Tech != "N/R")

# Recoding the answers & adding new column - here we we create a new variable that corresponds to the legislator preference on an issue-area.  
tech <- tech %>%
  mutate(Tech_num = recode(Tech,
                             "Em certas circunstâncias,  o STF tem o poder de decidir sobre temas de comportamento (aborto e casamento) e de revogar atos de ÓRGÃOS TÉCNICOS (agências reguladoras e instituições de controle)." = "Active Court",
                             "O STF NÃO tem o poder de decidir sobre temas de comportamento (aborto e casamento) e NEM de revogar atos ÓRGÃOS TÉCNICOS." = "Passive Court"))

# Calculating the Confidence Interval 

## Creating a count and percentage table
tabela_cenarios = tech %>% group_by(Tech_num) %>% summarise(count = n())
perct = tabela_cenarios$count/nrow(tech)
tabela_cenarios = cbind(tabela_cenarios, perct)


c1 = tabela_cenarios[1,3]
c2 = tabela_cenarios[2,3]

# Confidence intervals - here we calculate lower/upper confidence intervals
c1+c2
se = sqrt((c1*c2)/nrow(tech))
CI = 1.96 * se

lower.CI = c(c1-CI, c2-CI) 
upper.CI = c(c1+CI, c2+CI) 

# Confidence interval table - here we build the confidence interval table
tabela_cenarios_tech = cbind(tabela_cenarios, lower.CI, upper.CI)
tabela_cenarios_tech

# Printing the table
stargazer(tabela_cenarios_tech,summary = F, type = "text", digits=4, no.space=T, flip=F)

# Saving confidence interval table [Bureaucracy]
writexl::write_xlsx(tabela_cenarios_tech, "Bureaucracy.xlsx")

#### Figure 5: Bureaucracy

a = read_xlsx("Bureaucracy.xlsx")

graph_tech = ggplot(a) +
  geom_bar( aes(Tech_num , y=perct), stat="identity", fill="white", width=0.4, colour = "black", alpha= 2.5) +
  geom_errorbar( aes(x=Tech_num, ymin=lower.CI, ymax=upper.CI), width=0.1, colour="black", alpha=0.9, size=1) + ggtitle("Figure 5: Bureaucracy") +
  xlab("") +
  ylab("") +
  theme(plot.title = element_text(family = "serif", size = 20, face = "bold", hjust = 0.5), axis.title.x = element_text(family = "serif"), axis.text = element_text(family = "serif", colour="black", size = 19) ) 

ggsave("Figure 5 - Bureaucracy.png", width = 2105 , height = 2105, units="px", dpi= 420)

```

```{r}
#### Merging Images

# Create the combined plot
combined_plot <- ggarrange(
  graph_pork, 
  ggarrange(graph_check, graph_elec, graph_venue, graph_tech, ncol = 2, nrow = 2)
)

# Save the plot
ggsave("Figure 1 to 5.png", plot = combined_plot, width = 8500, height = 3500,  units = "px", dpi = 500)

```

```{r}
#### Figure 6: Power Analysis Over Sample Sizes (Earmark Distribution)

# Set a range of sample sizes
sample_sizes <- seq(20, 200, 10)

power_values <- sapply(sample_sizes, function(n) {
  result <- pwrss.z.prop(p = 0.3223, p0 = 0.5, alpha = 0.05, n = n, alternative = "not equal")
  return(result$power)
})

# Create a data frame for the results
data <- data.frame(N = sample_sizes, Statistical_Power = power_values)

# Plot the results using ggplot2
ggplot(data, aes(x = N, y = Statistical_Power)) +
  geom_line(color = "black") +
  geom_point(size = 1) +
  labs(title = "Earmark Distribution",
       x = "Sample Size (N)",
       y = "Statistical Power") +
  theme_minimal()

ggsave("Figure 6 - Power Analysis Over Sample Sizes (Earmark Distribution).png", units="px", dpi= 420)

```

```{r}
#### Figure 7: Distributions of z-statistic for Proportions (Earmark)

pwrss.graph = pwrss.z.prop(p = 0.6771, p0 = 0.5,
             alpha = 0.05, n = 96, 
             alternative = "not equal")

plot(pwrss.graph)  ## The image was saved manually as "Figure 7 - Distributions of z-statistic for Proportions (Earmark).png"
```


```{r}
# Table A.5: Legislative preferences by issue-area and scenario


# This code creates a summary data frame (`final_table`) for reporting proportions, confidence intervals, and sample sizes corresponding to five distinct issue areas (e.g., earmark distribution, social values). 

final_table <- data.frame(
  Figure = c("Figure 1", "", "Figure 2", "", "Figure 3", "", "Figure 4", "", "Figure 5", ""),
  `Issue-area` = c(
    "Earmark distribution", "", 
    "Social values", "", 
    "Electoral Legislation", "", 
    "Privileged jurisdiction", "", 
    "Bureaucracy", ""
  ),
  Scenario = rep(c("Active Court", "Passive Court"), 5),
  Proportion = c(0.3229, 0.6771, 0.4583, 0.5417, 0.4700, 0.5300, 0.4433, 0.5567, 0.4065, 0.5934),
  SE = c(0.0475, 0.0475, 0.0510, 0.0510, 0.0509, 0.0509, 0.0521, 0.0521, 0.0505, 0.0505),
  `95% CI (Lower)` = c("[0.2294, 0.4165]",
"[0.5835, 0.7706]",
"[0.3587, 0.5580]",
"[0.4420, 0.6413]",
"[0.3702, 0.5698]",
"[0.4302, 0.6298]",
"[0.3412, 0.5454]",
"[0.4546, 0.6588]",
"[0.3057, 0.5075]",
"[0.4925, 0.6943]"),
  `Sample Size (N)` = c(31, 65, 44, 52, 47, 53, 43, 54, 37, 54)
)


# Saving the table as HTML 
stargazer(
  final_table,
  type = "html",                          
  summary = FALSE, 
  rownames = FALSE,
  covariate.labels = c("Figure",	"Issue-area",	"Scenario",	"Proportion",	"SE",	"95% CI","N"),                      
  out = "Table A.5.html",            
  digit.separate = 0,                     
  digits = 4                              
)
```

```{r}
#### Table A.6: Power Analysis Summary

# Extracting the calculated values
power <- round(pwrss.graph$power, 3)               # Statistical Power (1 - β)
beta <- round(1 - pwrss.graph$power, 3)            # Type II Error Rate (β)
ncp <- round(pwrss.graph$ncp, 3)                   # Non-centrality Parameter
alpha <- pwrss.graph$parms$alpha                   # Type I Error Rate (α)
n <- pwrss.graph$n                                 # Sample Size (n)
alternative <- ifelse(pwrss.graph$parms$alternative == "not equal", 
                      "Not Equal", 
                      pwrss.graph$parms$alternative) # Alternative Hypothesis (formatted)

# Create the table as a data frame
parameter_table <- data.frame(
  Parameter = c(
    "Approach",                        # Row 1
    "Test Type",                       # Row 2
    "Null Hypothesis (H0)",            # Row 3
    "Alternative Hypothesis (HA)",     # Row 4
    "Statistical Power",               # Row 5
    "Sample Size (n)",                 # Row 6
    "Alternative Hypothesis",          # Row 7
    "Non-centrality Parameter",        # Row 8
    "Type I Error Rate (α)",           # Row 9
    "Type II Error Rate (β)"           # Row 10
  ),
  Value = c(
    "Normal Approximation",            # Row 1
    "Proportion against a Constant (z Test)", # Row 2
    "p = p0",                          # Row 3
    "p ≠ p0",                          # Row 4
    power,                             # Row 5: Dynamically calculated
    n,                                 # Row 6: Dynamically calculated
    alternative,                       # Row 7: Dynamically calculated
    ncp,                               # Row 8: Dynamically calculated
    alpha,                             # Row 9: Dynamically calculated
    beta                               # Row 10: Dynamically calculated
  ),
  stringsAsFactors = FALSE 
)

# Print the table to the console
print(parameter_table, row.names = FALSE)

# Better visual output
kable(parameter_table, format = "html", table.attr = 'class="table table-bordered"')


```




